import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'SimSun']
plt.rcParams['axes.unicode_minus'] = False

# 读取数据
file_path = r"D:\C题\处理后数据(1).xlsx"
try:
    data = pd.read_excel(file_path)
    print("数据读取成功！")
    print(f"数据形状: {data.shape}")
    print("\n前5行数据:")
    print(data.head())
    print("\n列名和数据类型:")
    print(data.dtypes)
except Exception as e:
    print(f"读取文件时出错: {e}")
    # 创建示例数据用于演示
    print("创建示例数据进行演示...")
    np.random.seed(42)
    n_samples = 100
    data = pd.DataFrame({
        '达标概率': np.random.normal(0.8, 0.15, n_samples),
        '潜在风险': np.random.normal(0.3, 0.1, n_samples),
        '关键变量1': np.random.normal(10, 2, n_samples),
        '关键变量2': np.random.normal(5, 1, n_samples)
    })

# 自动识别数值列的函数
def identify_numeric_columns(data):
    """识别数据中的数值列"""
    numeric_columns = []
    for col in data.columns:
        try:
            # 尝试转换为数值类型
            numeric_data = pd.to_numeric(data[col], errors='coerce')
            # 如果超过80%的数据可以转换为数值，则认为这是数值列
            if numeric_data.notna().mean() > 0.8:
                numeric_columns.append(col)
                print(f"识别为数值列: {col}")
        except:
            continue
    return numeric_columns

# 定义蒙特卡洛模拟函数
def monte_carlo_simulation(data, n_simulations=1000, error_std=0.1):
    """
    执行蒙特卡洛模拟，分析误差对达标概率和潜在风险的影响
    """
    results = {
        '达标概率': [],
        '潜在风险': [],
        '达标概率_mean': [],
        '潜在风险_mean': [],
        '达标概率_std': [],
        '潜在风险_std': []
    }
    
    # 自动识别数值列
    numeric_columns = identify_numeric_columns(data)
    
    if len(numeric_columns) < 2:
        print("错误: 数据中需要至少2个数值列")
        print("请手动指定数值列名称")
        return results
    
    # 使用前两个数值列
    target_col = numeric_columns[0]
    risk_col = numeric_columns[1]
    
    print(f"使用 '{target_col}' 作为达标指标")
    print(f"使用 '{risk_col}' 作为风险指标")
    
    # 确保数据是数值类型
    data[target_col] = pd.to_numeric(data[target_col], errors='coerce')
    data[risk_col] = pd.to_numeric(data[risk_col], errors='coerce')
    
    # 删除NaN值
    data_clean = data[[target_col, risk_col]].dropna()
    
    if len(data_clean) == 0:
        print("错误: 清理后没有有效数据")
        return results
    
    # 基准值（无误差）
    baseline_target = data_clean[target_col].mean()
    baseline_risk = data_clean[risk_col].mean()
    
    print(f"基准达标概率: {baseline_target:.4f}")
    print(f"基准潜在风险: {baseline_risk:.4f}")
    
    # 蒙特卡洛模拟
    for i in range(n_simulations):
        # 添加随机误差
        error = np.random.normal(0, error_std, len(data_clean))
        
        # 计算带误差的指标
        target_with_error = data_clean[target_col] + error
        risk_with_error = data_clean[risk_col] + error * 0.5
        
        # 计算统计量
        target_mean = target_with_error.mean()
        risk_mean = risk_with_error.mean()
        
        # 存储结果
        results['达标概率'].append(target_mean)
        results['潜在风险'].append(risk_mean)
        
        # 每100次模拟计算一次均值和标准差
        if (i + 1) % 100 == 0:
            results['达标概率_mean'].append(np.mean(results['达标概率'][-100:]))
            results['潜在风险_mean'].append(np.mean(results['潜在风险'][-100:]))
            results['达标概率_std'].append(np.std(results['达标概率'][-100:]))
            results['潜在风险_std'].append(np.std(results['潜在风险'][-100:]))
    
    # 添加基准值
    results['baseline_target'] = baseline_target
    results['baseline_risk'] = baseline_risk
    
    return results

# 执行蒙特卡洛模拟
print("\n开始蒙特卡洛模拟...")
simulation_results = monte_carlo_simulation(data, n_simulations=1000, error_std=0.05)

# 检查是否有有效结果
if not simulation_results['达标概率']:
    print("模拟失败，请检查数据格式")
else:
    # 创建图表
    fig, axes = plt.subplots(2, 2, figsize=(16, 12))

    # 1. 达标概率分布
    axes[0, 0].hist(simulation_results['达标概率'], bins=30, alpha=0.7, color='skyblue', edgecolor='black')
    axes[0, 0].axvline(simulation_results['baseline_target'], color='red', linestyle='--', linewidth=2, 
                      label=f'基准值: {simulation_results["baseline_target"]:.3f}')
    axes[0, 0].set_xlabel('达标概率', fontsize=12, fontweight='bold')
    axes[0, 0].set_ylabel('频数', fontsize=12, fontweight='bold')
    axes[0, 0].set_title('达标概率的蒙特卡洛模拟分布', fontsize=14, fontweight='bold', pad=20)
    axes[0, 0].legend(fontsize=10)
    axes[0, 0].grid(True, alpha=0.3)

    # 2. 潜在风险分布
    axes[0, 1].hist(simulation_results['潜在风险'], bins=30, alpha=0.7, color='lightcoral', edgecolor='black')
    axes[0, 1].axvline(simulation_results['baseline_risk'], color='red', linestyle='--', linewidth=2, 
                      label=f'基准值: {simulation_results["baseline_risk"]:.3f}')
    axes[0, 1].set_xlabel('潜在风险', fontsize=12, fontweight='bold')
    axes[0, 1].set_ylabel('频数', fontsize=12, fontweight='bold')
    axes[0, 1].set_title('潜在风险的蒙特卡洛模拟分布', fontsize=14, fontweight='bold', pad=20)
    axes[0, 1].legend(fontsize=10)
    axes[0, 1].grid(True, alpha=0.3)

    # 3. 达标概率随时间变化
    if simulation_results['达标概率_mean']:
        x_points = range(100, 1001, 100)
        axes[1, 0].plot(x_points, simulation_results['达标概率_mean'], 'o-', color='blue', label='均值', markersize=6)
        axes[1, 0].fill_between(x_points, 
                               np.array(simulation_results['达标概率_mean']) - np.array(simulation_results['达标概率_std']),
                               np.array(simulation_results['达标概率_mean']) + np.array(simulation_results['达标概率_std']),
                               alpha=0.2, color='blue', label='±1标准差')
        axes[1, 0].axhline(simulation_results['baseline_target'], color='red', linestyle='--', label='基准值')
        axes[1, 0].set_xlabel('模拟次数', fontsize=12, fontweight='bold')
        axes[1, 0].set_ylabel('达标概率', fontsize=12, fontweight='bold')
        axes[1, 0].set_title('达标概率随模拟次数的变化', fontsize=14, fontweight='bold', pad=20)
        axes[1, 0].legend(fontsize=10)
        axes[1, 0].grid(True, alpha=0.3)
        axes[1, 0].tick_params(axis='both', which='major', labelsize=10)

    # 4. 潜在风险随时间变化
    if simulation_results['潜在风险_mean']:
        x_points = range(100, 1001, 100)
        axes[1, 1].plot(x_points, simulation_results['潜在风险_mean'], 'o-', color='red', label='均值', markersize=6)
        axes[1, 1].fill_between(x_points, 
                               np.array(simulation_results['潜在风险_mean']) - np.array(simulation_results['潜在风险_std']),
                               np.array(simulation_results['潜在风险_mean']) + np.array(simulation_results['潜在风险_std']),
                               alpha=0.2, color='red', label='±1标准差')
        axes[1, 1].axhline(simulation_results['baseline_risk'], color='green', linestyle='--', label='基准值')
        axes[1, 1].set_xlabel('模拟次数', fontsize=12, fontweight='bold')
        axes[1, 1].set_ylabel('潜在风险', fontsize=12, fontweight='bold')
        axes[1, 1].set_title('潜在风险随模拟次数的变化', fontsize=14, fontweight='bold', pad=20)
        axes[1, 1].legend(fontsize=10)
        axes[1, 1].grid(True, alpha=0.3)
        axes[1, 1].tick_params(axis='both', which='major', labelsize=10)

    # 确保布局调整
    plt.tight_layout(pad=3.0)
    plt.show()

    # 输出统计结果
    print("\n=== 蒙特卡洛模拟结果统计 ===")
    print(f"模拟次数: 1000")
    print(f"达标概率 - 基准值: {simulation_results['baseline_target']:.4f}")
    print(f"达标概率 - 模拟均值: {np.mean(simulation_results['达标概率']):.4f}")
    print(f"达标概率 - 模拟标准差: {np.std(simulation_results['达标概率']):.4f}")
    print(f"达标概率 - 95%置信区间: [{np.percentile(simulation_results['达标概率'], 2.5):.4f}, "
          f"{np.percentile(simulation_results['达标概率'], 97.5):.4f}]")

    print(f"\n潜在风险 - 基准值: {simulation_results['baseline_risk']:.4f}")
    print(f"潜在风险 - 模拟均值: {np.mean(simulation_results['潜在风险']):.4f}")
    print(f"潜在风险 - 模拟标准差: {np.std(simulation_results['潜在风险']):.4f}")
    print(f"潜在风险 - 95%置信区间: [{np.percentile(simulation_results['潜在风险'], 2.5):.4f}, "
          f"{np.percentile(simulation_results['潜在风险'], 97.5):.4f}]")

    # 计算误差影响程度
    target_impact = (np.std(simulation_results['达标概率']) / simulation_results['baseline_target']) * 100
    risk_impact = (np.std(simulation_results['潜在风险']) / simulation_results['baseline_risk']) * 100

    print(f"\n误差影响程度:")
    print(f"达标概率的相对波动: {target_impact:.2f}%")
    print(f"潜在风险的相对波动: {risk_impact:.2f}%")

print("\n分析完成！")
